library(magrittr)
library(tidyverse)
library(Seurat)
library(readxl)
library(cowplot)
library(colorblindr)
library(viridis)
library(progeny)

theme_cowplot2 <- function(...) {
  theme_cowplot(font_size = 12, ...) %+replace%
    theme(strip.background = element_blank(),
          plot.background = element_blank())
}
theme_set(theme_cowplot2())
coi <- params$cell_type
cell_sort <- params$cell_sort
cell_type_major <- params$cell_type_major
louvain_resolution <- params$louvain_resolution
louvain_cluster <- params$louvain_cluster

1 Cluster markers

1.1 Major Ovarian.cancer.cell markers for cell assign

### load all data ---------------------------------
source("_src/global_vars.R")

# seu_obj <- read_rds(paste0("/work/shah/isabl_data_lake/analyses/16/52/1652/celltypes/", coi, "_processed.rds"))
seu_obj <- read_rds(paste0("/work/shah/uhlitzf/data/SPECTRUM/freeze/v6/Ovarian.cancer.cell_highqc.rds"))

myfeatures <- c("umapharmony_1", "umapharmony_2", "sample", louvain_cluster, "doublet", "nCount_RNA", "nFeature_RNA", "percent.mt", "doublet_score")

plot_data_wrapper <- function(cluster_res) {
  cluster_res <- enquo(cluster_res)
  as_tibble(FetchData(seu_obj, myfeatures)) %>% 
    left_join(meta_tbl, by = "sample") %>% 
    rename(cluster = !!cluster_res) %>% 
    mutate(cluster = as.character(cluster),
           tumor_supersite = ordered(tumor_supersite, levels = names(clrs$tumor_supersite)))
}

plot_data <- plot_data_wrapper(louvain_cluster)

WFDC2, CD24, CLDN3, KRT7, KRT8, KRT17, KRT18, KRT19, EPCAM, WT1, CLDN4, MSLN, FOLR1, MUC1

1.2 Subtype currated markers

helper_f <- function(x) ifelse(is.na(x), "", x)

markers_v6_super[[coi]] %>% 
  group_by(subtype) %>% 
  mutate(rank = row_number(gene)) %>% 
  spread(subtype, gene) %>% 
  mutate_all(.funs = helper_f) %>% 
  formattable::formattable()
rank Cancer.cell.1 Cancer.cell.2 Cancer.cell.3 Cancer.cell.4 Cancer.cell.5 Cancer.cell.6 Ciliated.cell.1 Ciliated.cell.2 Ciliated.cell.3 Cycling.cancer.cell.1 Cycling.cancer.cell.2 doublet.Immune.cell
1 DAPL1 IGFBP3 CXCL10 COL1A1 CCND1 FTL C20orf85 C5orf49 ATP5F1E CENPF ATAD2 CD52
2 FOLR1 KRT17 IFI6 COL1A2 CRISP3 HSPA6 CAPS CFAP126 DEFB1 HMGB2 H2AFZ IGHG1
3 SCGB1D2 KRT7 IFIT1 COL3A1 OVGP1 MGP CAPSL DNALI1 GNAS MKI67 HIST1H4C IGKC
4 SCGB2A1 LCN2 IFIT2 DCN PLCG2 RACK1 PIFO PIFO GPX3 PTTG1 PCLAF IGLC3
5 SNHG19 MMP7 IFIT3 FN1 PTMS SCX TPPP3 RSPH1 IMPA2 TOP2A TUBA1B IL7R
6 S100A10 ISG15 SPARC SNHG9 ZFAS1 PTPRC
7 S100A9 MX1
8 SLPI
9 TACSTD2

1.3 Subtype cluster markers

marker_tbl <- read_tsv("/work/shah/uhlitzf/data/SPECTRUM/freeze/v6/Ovarian.cancer.cell_highqc_markers_02.tsv")

## Hypergeometric test --------------------------------------

test_set <- marker_tbl %>% 
  group_by(cluster) %>% 
  filter(!str_detect(gene, "^RPS|^RPL")) %>% 
  slice(1:50) %>% 
  mutate(k = length(cluster)) %>% 
  ungroup %>%
  select(cluster, gene, k) %>% 
  mutate(join_helper = 1) %>% 
  group_by(cluster, join_helper, k) %>%
  nest(test_set = gene)

markers_doub_tbl <- markers_v6 %>% 
  enframe("subtype", "gene") %>% 
  filter(!(subtype %in% unique(c(coi, cell_type_major)))) %>% 
  unnest(gene) %>% 
  group_by(gene) %>% 
  filter(length(gene) == 1) %>% 
  mutate(subtype = paste0("doublet.", subtype)) %>% 
  bind_rows(tibble(subtype = "Mito.high", gene = grep("^MT-", rownames(seu_obj), value = T)))

ref_set <- markers_v6_super[[coi]] %>% 
  bind_rows(markers_doub_tbl) %>% 
  group_by(subtype) %>% 
  mutate(m = length(gene),
         n = length(rownames(seu_obj))-m,
         join_helper = 1) %>% 
  group_by(subtype, m, n, join_helper) %>%
  nest(ref_set = gene)

hyper_tbl <- test_set %>% 
  left_join(ref_set, by = "join_helper") %>% 
  group_by(cluster, subtype, m, n, k) %>%
  do(q = length(intersect(unlist(.$ref_set), unlist(.$test_set)))) %>%
  mutate(pval = 1-phyper(q = q, m = m, n = n, k = k)) %>%
  ungroup %>%
  mutate(qval = p.adjust(pval, "BH"),
         sig = qval < 0.01)

# hyper_tbl %>% 
#   group_by(subtype) %>% 
#   filter(any(qval < 0.01)) %>%
#   ggplot(aes(subtype, -log10(qval), fill = sig)) +
#   geom_bar(stat = "identity") +
#   facet_wrap(~cluster) +
#   coord_flip()
  
low_rank <- str_detect(unique(hyper_tbl$subtype), "doublet|dissociated")
subtype_lvl <- c(sort(unique(hyper_tbl$subtype)[!low_rank]), sort(unique(hyper_tbl$subtype)[low_rank]))
  
cluster_label_tbl <- hyper_tbl %>% 
  mutate(subtype = ordered(subtype, levels = subtype_lvl)) %>% 
  arrange(qval, subtype) %>%
  group_by(cluster) %>% 
  slice(1) %>% 
  mutate(subtype = ifelse(sig, as.character(subtype), paste0("unknown_", cluster))) %>% 
  select(cluster, cluster_label = subtype) %>% 
  ungroup %>% 
  mutate(cluster_label = make.unique(cluster_label, sep = "_"))

seu_obj$cluster_label <- unname(deframe(cluster_label_tbl)[as.character(unlist(seu_obj[[paste0("RNA_snn_res.", louvain_resolution)]]))])
plot_data$cluster_label <- seu_obj$cluster_label

cluster_n_tbl <- seu_obj$cluster_label %>% 
  table() %>% 
  enframe("cluster_label", "cluster_n") %>% 
  mutate(cluster_nrel = cluster_n/sum(cluster_n))

marker_sheet <- marker_tbl %>% 
  left_join(cluster_label_tbl, by = "cluster") %>% 
  group_by(cluster_label) %>% 
  filter(!str_detect(gene, "^RPS|^RPL")) %>% 
  slice(1:50) %>% 
  mutate(rank = row_number(-avg_logFC)) %>% 
  select(cluster_label, gene, rank) %>% 
  ungroup %>% 
  mutate(cluster_label = ordered(cluster_label, levels = unique(c(names(clrs$cluster_label[[coi]]), sort(cluster_label))))) %>% 
  spread(cluster_label, gene) %>% 
  mutate_all(.funs = helper_f)

marker_tbl_annotated <- marker_tbl %>% 
  left_join(cluster_label_tbl, by = "cluster") %>% 
  left_join(cluster_n_tbl, by = "cluster_label") %>% 
  select(-cluster) %>% 
  mutate(cluster_label = ordered(cluster_label, levels = unique(c(names(clrs$cluster_label[[coi]]), sort(cluster_label))))) %>% 
  arrange(cluster_label, -avg_logFC, p_val_adj)

write_tsv(marker_sheet, paste0("/work/shah/uhlitzf/data/SPECTRUM/freeze/v6/", coi, "_marker_sheet.tsv"))
write_tsv(marker_sheet, paste0("/work/shah/uhlitzf/data/SPECTRUM/freeze/v6/supplementary_tables/", coi, "_marker_sheet.tsv"))

write_tsv(marker_tbl_annotated, paste0("/work/shah/uhlitzf/data/SPECTRUM/freeze/v6/", coi, "_marker_table_annotated.tsv"))
write_tsv(marker_tbl_annotated, paste0("/work/shah/uhlitzf/data/SPECTRUM/freeze/v6/supplementary_tables/", coi, "_marker_table_annotated.tsv"))

formattable::formattable(marker_sheet)
rank Cancer.cell.1 Cancer.cell.2 Cancer.cell.3 Cancer.cell.4 Cancer.cell.5 Cancer.cell.6 Cycling.cancer.cell.1 Cycling.cancer.cell.2 Ciliated.cell.1 Ciliated.cell.2 Ciliated.cell.3 doublet.Immune.cell
1 DAPL1 S100A9 ISG15 COL1A1 OVGP1 HSPA6 CENPF HIST1H4C CAPS CFAP126 GPX3 IGKC
2 SCGB2A1 KRT17 CXCL10 COL3A1 PLCG2 FTL PTTG1 TUBA1B C20orf85 PIFO GNAS IGLC3
3 SCGB1D2 LCN2 IFIT3 COL1A2 CRISP3 RACK1 UBE2S PCLAF TPPP3 DNALI1 ATP5F1E LTB
4 SNHG19 TACSTD2 IFIT1 FN1 SNHG9 ZFAS1 HMGB2 ATAD2 C9orf24 RSPH1 DEFB1 IGHG1
5 FOLR1 IGFBP3 MX1 SPARC CCND1 SCX TPX2 H2AFZ RSPH1 C5orf49 IMPA2 IL7R
6 CRABP1 SLPI IFI6 DCN PTMS MGP ASPM PCNA C5orf49 LRRIQ1 GSTM3 CD52
7 PTGDS S100A10 IFIT2 LUM PEG10 IGFBP7 MKI67 CKS1B TMEM190 FAM183A PSMA7 PTPRC
8 THSD4 KRT7 IFI44L IGFBP7 DSTN NPM1 CCNB1 TK1 PIFO C9orf24 ANXA4 GIMAP7
9 APOA1 RARRES1 PARP14 COL6A1 YWHAE CNBP CDC20 MKI67 FAM183A CASC1 FTL SRGN
10 GPRC5A TNFSF10 COL6A2 SNHG19 FKBP4 TOP2A TYMS C1orf194 C20orf85 CAPS IGHG3
11 MMP7 STAT1 CCDC80 APOLD1 EEF1D CKS2 TUBB CETN2 ARMC3 IGFBP7 IGLC2
12 CRYAB IFI27 CTHRC1 MYL12B CDV3 UBE2C HMGB2 ODF3B C1orf194 PCLAF BTG1
13 FTH1 RSAD2 CALD1 PTPRA GAPDH ARL6IP1 DUT TXN TPPP3 PLIN2 CCR7
14 NDRG1 OAS1 VIM WBP11 OAZ1 BIRC5 TOP2A AGR3 CFAP45 PFDN4 KLF2
15 KRT19 XAF1 COL6A3 SNHG25 YBX3 LGALS1 DEK CFAP126 EFCAB1 C12orf75 BIRC3
16 S100A11 GBP1 NNMT RBMX YBX1 CCNB2 SMC4 IGFBP7 DNAAF1 RHOBTB3 CD3D
17 VEGFA SAMD9 SPARCL1 EPB41L2 LTC4S TUBA1B MCM7 CAPSL ENKUR CDKN2C RGS1
18 ANXA1 IFIH1 LGALS1 RBM8A CYC1 CEP55 UBE2C C11orf88 CEP126 BMP7 SPP1
19 FTL GBP4 MGP PABPC4 NOP53 CDKN3 HELLS MORN2 HYDIN EIF3K IL32
20 SLC2A1 PSMB9 TAGLN METAP2 EEF2 NUSAP1 CDK1 TUBB4B PPIL6 XAGE3 CD48
21 SOX4 DDX58 MMP11 AMD1 EDN1 JPT1 MCM3 TUBA1A CAPSL DSTN FYB1
22 S100A6 LY6E VCAN STRAP EIF4A2 H2AFZ CENPF LRRIQ1 ZMYND10 MSRB1 TRBC2
23 TFPI2 ISG20 ACTA2 BRD9 TSTA3 NMU STMN1 FOXJ1 TEKT2 EIF1 IKZF1
24 PI3 OAS3 RARRES2 CREBZF PRELID1 KPNA2 MYBL2 ZMYND10 CFAP43 ATP5PO CORO1A
25 C15orf48 HLA-B COL5A2 MYL6B SLU7 HMMR CLSPN DNAAF1 DYDC2 RDH10 CD3G
26 NEAT1 CXCL11 AEBP1 PPT1 CANX PTMS UBE2T C9orf116 CFAP44 COX6B1 B2M
27 S100A2 RNF213 POSTN SNRPB CTTNBP2 HMGN2 NASP EFCAB1 CDC20B SDC2 CD2
28 TNFRSF12A LAP3 PLA2G2A MYL12A TTC1 ECT2 MCM4 EFHC1 CFAP54 SFTA2 GMFG
29 C19orf33 IFI16 ERP27 TMCO1 CLINT1 SMC4 SMC2 DNALI1 FAM81B ORM1 RIPOR2
30 TCIM TAP1 NBL1 SNRPB2 TMEM238 HMGB1 NUSAP1 DYDC2 SPATA17 CNTD2 ETS1
31 TAGLN IDO1 COL8A1 NAP1L1 COX4I1 CENPE HMGB1 IK TEKT1 HOXB-AS3 GIMAP4
32 MAL2 HLA-A TIMP3 ST13 ZNF787 STMN1 RAD51AP1 CEP126 SPEF1 PTH1R HCST
33 LGALS3 PLSCR1 LRRC75A PTX3 CLTB ANP32E ZWINT SNTN CFAP46 CYS1 TRAC
34 NUPR1 HERC5 IGFBP4 MIR4458HG HIGD2A CKS1B TPX2 AL357093.2 DNAH7 GLRX TRBC1
35 ANXA2 IFI44 CDH6 WDR70 NACA KIF20B DNMT1 ROPN1L CFAP73 EEF1A2 TMSB4X
36 CST6 B2M MMP2 WDR45B UQCRB RAD21 RRM2 MLF1 ZBBX ERVMER34-1 CXCR4
37 EMP1 HLA-C BGN PA2G4 ATP5PB NUCKS1 E2F1 PSENEN ROPN1L CA2 LAPTM5
38 IL32 OAS2 TPM1 NMT1 RSPO3 HMGB3 CENPW ERICH3 EFCAB10 LINC02532 CD3E
39 SQSTM1 BST2 IGFBP6 TTYH1 RSL1D1 TROAP RANBP1 CCDC170 WDR38 NCCRP1 C1QB
40 CAST WARS C1R VIM CYHR1 DLGAP5 MAD2L1 TEKT2 SPAG8 LAMB4 LCP1
41 MUC4 SAMD9L SPON2 KRT81 UBE2D2 AURKA TMPO PERP DRC1 BICC1 ARHGDIB
42 TGM2 EIF2AK2 COL5A1 IGFBP3 S100A1 H2AFV CENPX SPA17 CFAP157 MIOX CD37
43 SAT1 UBE2L6 KRT7 NSA2 SRP72 TUBB4B FANCI C9orf135 RSPH4A ORM2 TYROBP
44 CLDN4 MX2 HTRA1 HSP90B1 SSB CENPA TMEM106C LRRC23 WDR49 CPNE8 CCL4
45 FAM3C OASL CAV1 HMGB1 PUF60 NCAPD2 DHFR CFAP45 C7orf57 NID2 EMP3
46 BIRC3 IFI35 SYT8 TMA7 BTF3 CENPW PRKDC ARMC3 DYNLRB2 AC022784.1 LSP1
47 LDHA IFITM3 C1S ANO1 UQCRH TUBA1C HMGN2 WDR54 WDR63 SLCO4A1 IGHG4
48 CAMK2N1 HLA-E SERPINF1 PEG3 EIF3E NEK2 BIRC5 MNS1 MAP3K19 CR1L TRAF3IP3
49 ITGB8 RARRES3 ANXA2 CCND2 MAP2K2 NUF2 DTYMK TEKT1 C6orf118 RBP7 STK17B
50 MUC16 NUPR1 CFH CD177 LINC02308 MAD2L1 ASPM TUBA4B DNAH12 PCAT19 ACAP1

2 Clusters

2.1 sizes

enframe(sort(table(seu_obj$cluster_label))) %>% 
  mutate(name = ordered(name, levels = rev(name))) %>% 
  ggplot() +
  geom_bar(aes(name, value), stat = "identity") +
  theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5)) +
  labs(y = c("#cells"), x = "cluster")

2.2 UMAP

alpha_lvl <- ifelse(nrow(plot_data) < 20000, 0.2, 0.1)
pt_size <- ifelse(nrow(plot_data) < 20000, 0.2, 0.05)

common_layers_disc <- list(  
  geom_point(size = pt_size, alpha = alpha_lvl),
  NoAxes(),
  guides(color = guide_legend(override.aes = list(size = 2, alpha = 1))),
  labs(color = "")
)

common_layers_cont <- list(  
  geom_point(size = pt_size, alpha = alpha_lvl),
  NoAxes(),
  scale_color_gradientn(colors = viridis(9)),
  guides(color = guide_colorbar())
)

ggplot(plot_data, aes(umapharmony_1, umapharmony_2, color = cluster_label)) + 
  common_layers_disc +
  #facet_wrap(~therapy) +
  ggtitle("Sub cluster")

3 Filtering out doublet and mito high clusters

3.1 UMAP

my_subtypes <- names(clrs$cluster_label[[coi]])

my_subtypes <- c(my_subtypes, unlist(lapply(paste0("_", 1:3), function(x) paste0(my_subtypes, x)))) %>% .[!str_detect(., "doublet")]

cells_to_keep <- colnames(seu_obj)[seu_obj$cluster_label %in% my_subtypes]
# seu_obj_sub <- subset(seu_obj, cells = cells_to_keep)
# seu_obj_sub <- RunUMAP(seu_obj_sub, dims = 1:20, reduction = "harmony", reduction.name = "umapharmony", reduction.key = "umapharmony_")
# seu_obj_sub$cluster_label <- seu_obj$cluster_label[colnames(seu_obj) %in% colnames(seu_obj_sub)]
# write_rds(seu_obj_sub, paste0("/work/shah/uhlitzf/data/SPECTRUM/freeze/v6/", coi, "_processed_filtered.rds"))
seu_obj_sub <- read_rds(paste0("/work/shah/uhlitzf/data/SPECTRUM/freeze/v6/", coi, "_processed_filtered.rds"))

plot_data_sub <- as_tibble(FetchData(seu_obj_sub, c(myfeatures, "cluster_label"))) %>% 
  left_join(meta_tbl, by = "sample") %>% 
  mutate(cell_id = colnames(seu_obj_sub),
         cluster_label = ordered(cluster_label, levels = my_subtypes))
  
if (cell_sort == "CD45+") {
plot_data_sub <- filter(plot_data_sub, sort_parameters != "singlet, live, CD45-", !is.na(tumor_supersite))
}

if (cell_sort == "CD45-") {
plot_data_sub <- filter(plot_data_sub, sort_parameters != "singlet, live, CD45+", !is.na(tumor_supersite))
}

ggplot(plot_data_sub, aes(umapharmony_1, umapharmony_2, color = cluster_label)) + 
  common_layers_disc +
  scale_color_manual(values = clrs$cluster_label[[coi]]) +
  #facet_wrap(~cluster_label) +
  ggtitle("Sub cluster") 

ggplot(plot_data_sub, aes(umapharmony_1, umapharmony_2, color = patient_id_short)) + 
  common_layers_disc +
  ggtitle("Patient") +
  # facet_wrap(~therapy) +
  scale_color_manual(values = clrs$patient_id_short)

ggplot(plot_data_sub, aes(umapharmony_1, umapharmony_2, color = tumor_supersite)) + 
  # geom_point(aes(umapharmony_1, umapharmony_2), 
  #            color = "grey90", size = 0.01, 
  #            data = select(plot_data_sub, -tumor_supersite)) +
  common_layers_disc +
  ggtitle("Site") +
  # facet_wrap(~therapy) +
  scale_color_manual(values = clrs$tumor_supersite)

write_tsv(select(plot_data_sub, cell_id, everything(), -contains("RNA_")), paste0("/work/shah/uhlitzf/data/SPECTRUM/freeze/v6/", coi, "_embedding.tsv"))

3.2 add module scores and pathway scores

signature_modules <- read_excel("_data/small/signatures/SPECTRUM v5 sub cluster markers.xlsx", sheet = 2, skip = 1, range = "M2:P100") %>% 
  gather(module, gene) %>% 
  na.omit() %>% 
  group_by(module) %>% 
  do(gene = c(.$gene)) %>% 
  {setNames(.$gene, .$module)}

signature_modules$ISG.module <- c("CCL5", "CXCL10", "IFNA1", "IFNB1", "ISG15", "IFI27L2", "SAMD9L")

# ## compute expression module scores
# for (i in 1:length(signature_modules)) {
#   seu_obj_sub <- AddModuleScore(seu_obj_sub, features = signature_modules[i], name = names(signature_modules)[i])
#   seu_obj_sub[[names(signature_modules)[i]]] <- seu_obj_sub[[paste0(names(signature_modules)[i], "1")]]
#   seu_obj_sub[[paste0(names(signature_modules)[i], "1")]] <- NULL
#   print(paste(names(signature_modules)[i], "DONE"))
# }
# 
# ## compute progeny scores
# progeny_list <- seu_obj_sub@assays$RNA@data[VariableFeatures(seu_obj_sub),] %>% 
#   as.matrix %>% 
#   progeny %>% 
#   as.data.frame %>% 
#   as.list
# 
# names(progeny_list) <- make.names(paste0(names(progeny_list), ".pathway"))
# 
# for (i in 1:length(progeny_list)) {
#   seu_obj_sub <- AddMetaData(seu_obj_sub, metadata = progeny_list[[i]], 
#                              col.name = names(progeny_list)[i])
# }
# 
# write_rds(seu_obj_sub, paste0("/work/shah/uhlitzf/data/SPECTRUM/freeze/v6/", coi, "_processed_filtered.rds"))
seu_obj_sub <- read_rds(paste0("/work/shah/uhlitzf/data/SPECTRUM/freeze/v6/", coi, "_processed_filtered.rds"))

3.3 marker heatmap

marker_top_tbl <- marker_sheet[,-1] %>% 
  slice(1:10) %>% 
  as.list %>% 
  .[!str_detect(names(.), "doublet")] %>% 
  enframe("cluster_label_x", "gene") %>% 
  unnest(gene)

plot_data_markers <- as_tibble(FetchData(seu_obj_sub, c("cluster_label", myfeatures, unique(marker_top_tbl$gene)))) %>% 
  gather(gene, value, -c(1:(length(myfeatures)+1))) %>% 
  left_join(meta_tbl, by = "sample") %>% 
  mutate(cluster_label = ordered(cluster_label, levels = my_subtypes)) %>% 
  group_by(cluster_label, gene) %>% 
  summarise(value = mean(value, na.rm = T)) %>% 
  group_by(gene) %>% 
  mutate(value = scales::rescale(value)) %>% 
  left_join(marker_top_tbl, by = "gene") %>% 
  mutate(cluster_label_x = ordered(cluster_label_x, levels = rev(names(clrs$cluster_label[[coi]]))))

ggplot(plot_data_markers) +
  geom_tile(aes(gene, cluster_label, fill = value)) +
  theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5)) +
  facet_grid(~cluster_label_x, scales = "free", space = "free") +
  scale_fill_gradientn(colors = viridis(9)) +
  labs(fill = "Scaled\nexpression") +
  theme(aspect.ratio = 1,
        axis.line = element_blank(),
        axis.ticks = element_blank(),
        axis.title = element_blank())

# ggsave(paste0("_fig/002_marker_heatmap_", coi, ".pdf"), width = nrow(marker_top_tbl)/6, height = 5)

3.4 composition

3.4.1 per site

comp_site_tbl <- plot_data_sub %>%
  filter(!is.na(tumor_supersite)) %>% 
  group_by(cluster_label, tumor_supersite) %>%
  tally %>%
  group_by(tumor_supersite) %>%
  mutate(nrel = n/sum(n)*100) %>%
  ungroup

pnrel_site <- ggplot(comp_site_tbl) +
  geom_bar(aes(tumor_supersite, nrel, fill = cluster_label),
           stat = "identity", position = position_stack()) +
  scale_y_continuous(expand = c(0, 0)) +
  theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5)) +
  scale_fill_manual(values = clrs$cluster_label[[coi]]) +
  labs(fill = "Cluster", y = "Fraction [%]", x = "")

pnabs_site <- ggplot(comp_site_tbl) +
  geom_bar(aes(tumor_supersite, n, fill = cluster_label),
           stat = "identity", position = position_stack()) +
  scale_y_continuous(expand = c(0, 0)) +
  theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5)) +
  scale_fill_manual(values = clrs$cluster_label[[coi]]) +
  labs(fill = "Cluster", y = "# cells", x = "") 

plot_grid(pnabs_site, pnrel_site, ncol = 2, align = "h")

# ggsave(paste0("_fig/02_deep_dive_", coi, "_comp_site.pdf"), width = 8, height = 4)

3.4.2 per sample

comp_tbl_sample_sort <- plot_data_sub %>% 
  group_by(tumor_subsite, tumor_supersite, patient_id, consensus_signature, therapy, cluster_label) %>% 
  tally %>% 
  group_by(tumor_subsite, tumor_supersite, patient_id, consensus_signature, therapy) %>% 
  mutate(nrel = n/sum(n)*100,
         nsum = sum(n),
         log10n = log10(n),
         label_supersite = "Site",
         label_mutsig = "Signature",
         label_therapy = "Rx") %>% 
  ungroup %>% 
  arrange(desc(therapy), tumor_supersite) %>% 
  mutate(tumor_subsite_rx = paste0(tumor_subsite, "_", therapy)) %>% 
  mutate(tumor_subsite = ordered(tumor_subsite, levels = unique(tumor_subsite)),
         tumor_subsite_rx = ordered(tumor_subsite_rx, levels = unique(tumor_subsite_rx))) %>% 
  arrange(patient_id) %>% 
  mutate(label_patient_id = ifelse(as.logical(as.numeric(fct_inorder(as.character(patient_id)))%%2), "Patient1", "Patient2"))

sample_id_x_tbl <- plot_data_sub %>% 
  mutate(sort_short_x = cell_sort) %>% 
  distinct(patient_id, sort_short_x, tumor_subsite, therapy, sample) %>% 
  unite(sample_id_x, patient_id, sort_short_x, tumor_subsite, therapy) %>% 
  arrange(sample_id_x)

comp_tbl_sample_sort %>% 
  mutate(sort_short_x = cell_sort) %>% 
  unite(sample_id_x, patient_id, sort_short_x, tumor_subsite_rx) %>% 
  select(sample_id_x, cluster_label, n, nrel, nsum) %>% 
  left_join(sample_id_x_tbl, by = "sample_id_x") %>% 
  write_tsv(paste0("/work/shah/uhlitzf/data/SPECTRUM/freeze/v6/", coi, "_subtype_compositions.tsv"))
ybreaks <- c(500, 1000, 2000, 4000, 6000, 8000, 10000, 15000, 20000)

max_cells_per_sample <- max(comp_tbl_sample_sort$nsum)
ymaxn <- ybreaks[max_cells_per_sample < ybreaks][1]

comp_plot_wrapper <- function(y = "nrel", switch = NULL) {
  if (y == "nrel") ylab <- paste0("Fraction\nof cells [%]")
  if (y == "n") ylab <- paste0("Number\nof cells")
  p <- ggplot(comp_tbl_sample_sort, 
              aes_string("tumor_subsite_rx", y, fill = "cluster_label")) + 
    facet_grid(~patient_id, space = "free", scales = "free", switch = switch) +
    coord_cartesian(clip = "off") + 
    scale_fill_manual(values = clrs$cluster_label[[coi]]) + 
    theme(axis.text.x = element_blank(),
          axis.title.y = element_text(angle = 0, vjust = 0.5, hjust = 0.5, 
                                      margin = margin(0, -0.4, 0, 0, unit = "npc")),
          axis.ticks.x = element_blank(),
          axis.title.x = element_blank(),
          axis.line.x = element_blank(),
          strip.text.y = element_blank(),
          strip.text.x = element_blank(),
          strip.background.y = element_blank(),
          strip.background.x = element_blank(),
          plot.margin = margin(0, 0, 0, 0, "npc")) + 
    labs(x = "",
         y = ylab) +
    guides(fill = FALSE)
  if (y == "nrel") p <- p + 
    geom_bar(stat = "identity") +
    scale_y_continuous(expand = c(0, 0), 
                       breaks = c(0, 50, 100), 
                       labels = c("0", "50", "100"))
  if (y == "n") p <- p + 
    geom_bar(stat = "identity", position = position_stack()) +
    scale_y_continuous(expand = c(0, 0), 
                       breaks = c(0, ymaxn/2, ymaxn),
                       limits = c(0, ymaxn),
                       labels = c("", ymaxn/2, ymaxn)) +
    expand_limits(y = c(0, ymaxn)) +
    theme(panel.grid.major.y = element_line(linetype = 1, color = "grey90", size = 0.5))
  return(p)
} 

common_label_layers <- list(
  geom_tile(color = "white", size = 0.15),
  theme(axis.text.x = element_blank(),
        axis.ticks = element_blank(),
        axis.title.x = element_blank(),
        axis.line.x = element_blank(),
        strip.text = element_blank(),
        strip.background = element_blank(),
        plot.margin = margin(0, 0, 0, 0, "npc")),
  scale_y_discrete(expand = c(0, 0)),
  labs(y = ""),
  guides(fill = FALSE),
  facet_grid(~patient_id, 
             space = "free", scales = "free")
)

comp_label_site <- ggplot(distinct(comp_tbl_sample_sort, tumor_subsite_rx, therapy, tumor_supersite, label_supersite, patient_id), 
       aes(tumor_subsite_rx, label_supersite, 
           fill = tumor_supersite)) + 
  scale_fill_manual(values = clrs$tumor_supersite) +
  common_label_layers

comp_label_rx <- ggplot(distinct(comp_tbl_sample_sort, tumor_subsite_rx, therapy, tumor_supersite, label_therapy, consensus_signature, patient_id), 
       aes(tumor_subsite_rx, label_therapy, 
           fill = therapy)) + 
  scale_fill_manual(values = c(`post-Rx` = "gold3", `pre-Rx` = "steelblue")) +
  common_label_layers

comp_label_mutsig <- ggplot(distinct(comp_tbl_sample_sort, tumor_subsite_rx, therapy, tumor_supersite, label_mutsig, consensus_signature, patient_id), 
       aes(tumor_subsite_rx, label_mutsig, 
           fill = consensus_signature)) + 
  scale_fill_manual(values = clrs$consensus_signature) +
  common_label_layers

patient_label_tbl <- distinct(comp_tbl_sample_sort, patient_id, .keep_all = T)

comp_label_patient_id <- ggplot(patient_label_tbl, aes(tumor_subsite_rx, label_patient_id)) + 
  scale_fill_manual(values = clrs$consensus_signature) +
  geom_text(aes(tumor_subsite_rx, label_patient_id, label = patient_id)) +
  facet_grid(~patient_id, 
             space = "free", scales = "free") +
  coord_cartesian(clip = "off") + 
  theme_void() +
  theme(strip.text = element_blank(),
        strip.background = element_blank(),
        plot.margin = margin(0, 0, 0, 0, "npc"))

hist_plot_wrapper <- function(x = "nrel") {
  if (x == "nrel") {
    xlab <- paste0("Fraction of cells [%]")
    bw <- 5
  }
  if (x == "log10n") {
    xlab <- paste0("Number of cells")
    bw <- 0.2
  }
  p <- ggplot(comp_tbl_sample_sort) +
    ggridges::geom_density_ridges(
      aes_string(x, "cluster_label", fill = "cluster_label"), color = "black",
      stat = "binline", binwidth = bw, scale = 3) +
    facet_grid(label_supersite~., 
               space = "free", scales = "free") +
    scale_fill_manual(values = clrs$cluster_label[[coi]]) + 
    theme(axis.text.y = element_blank(),
          axis.ticks.y = element_blank(),
          axis.title.y = element_blank(),
          axis.line.y = element_blank(),
          strip.text = element_blank(),
          strip.background = element_blank(),
          plot.margin = margin(0, 0, 0, 0, "npc")) +
    scale_x_continuous(expand = c(0, 0)) +
    scale_y_discrete(expand = c(0, 0)) +
    guides(fill = FALSE) +
    labs(x = xlab)
  if (x == "log10n") p <- p + expand_limits(x = c(0, 3)) + 
    scale_x_continuous(expand = c(0, 0), 
                       labels = function(x) parse(text = paste("10^", x)))
  return(p)
}

pcomp1 <- comp_plot_wrapper("n")
pcomp2 <- comp_plot_wrapper("nrel")

pcomp_grid <- plot_grid(comp_label_patient_id, 
                        pcomp1, pcomp2, 
                        comp_label_site, comp_label_rx, comp_label_mutsig,
                        ncol = 1, align = "v", 
                        rel_heights = c(0.15, 0.33, 0.33, 0.06, 0.06, 0.06))

phist1 <- hist_plot_wrapper("log10n")

pcomp_hist_grid <- ggdraw() +
  draw_plot(pcomp_grid, x = 0.01, y = 0, width = 0.85, height = 1) +
  draw_plot(phist1, x = 0.87, y = 0.05, width = 0.12, height = 0.8)

pcomp_hist_grid

# ggsave(paste0("_fig/02_composition_v6_",coi,".pdf"), pcomp_hist_grid, width = 10, height = 2)

3.4.3 site specific cluster enrichment

comp_tbl_z <- comp_tbl_sample_sort %>% 
  filter(therapy == "pre-Rx",
         !(tumor_supersite %in% c("Ascites", "Other"))) %>% 
  group_by(patient_id, cluster_label) %>% 
  arrange(patient_id, cluster_label, nrel) %>% 
  mutate(rank = row_number(nrel),
         z_rank = scales::rescale(rank)) %>% 
  mutate(mean_nrel = mean(nrel, na.rm = T),
         sd_nrel = sd(nrel, na.rm = T),
         z_nrel = (nrel - mean_nrel) / sd_nrel) %>% 
  ungroup()

ggplot(comp_tbl_z) +
  geom_boxplot(aes(tumor_supersite, z_nrel, color = tumor_supersite), 
               outlier.shape = NA) +
  geom_point(aes(tumor_supersite, z_nrel, color = tumor_supersite), position = position_jitter(), size = 0.1) +
  facet_grid(~cluster_label, scales = "free_x", space = "free_x") +
  theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5),
        aspect.ratio = 5) +
  scale_color_manual(values = clrs$tumor_supersite)

ggplot(comp_tbl_z) +
  geom_boxplot(aes(tumor_supersite, z_rank, color = tumor_supersite), 
               outlier.shape = NA) +
  geom_point(aes(tumor_supersite, z_rank, color = tumor_supersite), position = position_jitter(), size = 0.1) +
  facet_grid(~cluster_label, scales = "free_x", space = "free_x") +
  theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5),
        aspect.ratio = 5) +
  scale_color_manual(values = clrs$tumor_supersite)

4 session info

devtools::session_info()
## ─ Session info ───────────────────────────────────────────────────────────────
##  setting  value                       
##  version  R version 3.6.2 (2019-12-12)
##  os       Debian GNU/Linux 10 (buster)
##  system   x86_64, linux-gnu           
##  ui       X11                         
##  language (EN)                        
##  collate  en_US.UTF-8                 
##  ctype    en_US.UTF-8                 
##  tz       Etc/UTC                     
##  date     2021-01-20                  
## 
## ─ Packages ───────────────────────────────────────────────────────────────────
##  package      * version    date       lib
##  ape            5.3        2019-03-17 [2]
##  assertthat     0.2.1      2019-03-21 [2]
##  backports      1.1.10     2020-09-15 [1]
##  bibtex         0.4.2.2    2020-01-02 [2]
##  Biobase        2.46.0     2019-10-29 [2]
##  BiocGenerics   0.32.0     2019-10-29 [2]
##  bitops         1.0-6      2013-08-17 [2]
##  broom          0.7.2      2020-10-20 [1]
##  callr          3.4.2      2020-02-12 [1]
##  caTools        1.17.1.4   2020-01-13 [2]
##  cellranger     1.1.0      2016-07-27 [2]
##  cli            2.0.2      2020-02-28 [1]
##  cluster        2.1.0      2019-06-19 [3]
##  codetools      0.2-16     2018-12-24 [3]
##  colorblindr  * 0.1.0      2020-01-13 [2]
##  colorspace   * 1.4-2      2019-12-29 [2]
##  cowplot      * 1.0.0      2019-07-11 [2]
##  crayon         1.3.4      2017-09-16 [1]
##  data.table     1.12.8     2019-12-09 [2]
##  DBI            1.1.0      2019-12-15 [2]
##  dbplyr         2.0.0      2020-11-03 [1]
##  desc           1.2.0      2018-05-01 [2]
##  devtools       2.2.1      2019-09-24 [2]
##  digest         0.6.25     2020-02-23 [1]
##  dplyr        * 1.0.2      2020-08-18 [1]
##  ellipsis       0.3.1      2020-05-15 [1]
##  evaluate       0.14       2019-05-28 [2]
##  fansi          0.4.1      2020-01-08 [2]
##  farver         2.0.3      2020-01-16 [1]
##  fitdistrplus   1.0-14     2019-01-23 [2]
##  forcats      * 0.5.0      2020-03-01 [1]
##  formattable    0.2.0.1    2016-08-05 [1]
##  fs             1.5.0      2020-07-31 [1]
##  future         1.15.1     2019-11-25 [2]
##  future.apply   1.4.0      2020-01-07 [2]
##  gbRd           0.4-11     2012-10-01 [2]
##  gdata          2.18.0     2017-06-06 [2]
##  generics       0.0.2      2018-11-29 [2]
##  ggplot2      * 3.3.2      2020-06-19 [1]
##  ggrepel        0.8.1      2019-05-07 [2]
##  ggridges       0.5.2      2020-01-12 [2]
##  globals        0.12.5     2019-12-07 [2]
##  glue           1.3.2      2020-03-12 [1]
##  gplots         3.0.1.2    2020-01-11 [2]
##  gridExtra      2.3        2017-09-09 [2]
##  gtable         0.3.0      2019-03-25 [2]
##  gtools         3.8.1      2018-06-26 [2]
##  haven          2.3.1      2020-06-01 [1]
##  hms            0.5.3      2020-01-08 [1]
##  htmltools      0.4.0      2019-10-04 [2]
##  htmlwidgets    1.5.1      2019-10-08 [2]
##  httr           1.4.2      2020-07-20 [1]
##  ica            1.0-2      2018-05-24 [2]
##  igraph         1.2.5      2020-03-19 [1]
##  irlba          2.3.3      2019-02-05 [2]
##  jsonlite       1.7.1      2020-09-07 [1]
##  KernSmooth     2.23-16    2019-10-15 [3]
##  knitr          1.26       2019-11-12 [2]
##  labeling       0.3        2014-08-23 [2]
##  lattice        0.20-38    2018-11-04 [3]
##  lazyeval       0.2.2      2019-03-15 [2]
##  leiden         0.3.1      2019-07-23 [2]
##  lifecycle      0.2.0      2020-03-06 [1]
##  listenv        0.8.0      2019-12-05 [2]
##  lmtest         0.9-37     2019-04-30 [2]
##  lsei           1.2-0      2017-10-23 [2]
##  lubridate      1.7.9.2    2020-11-13 [1]
##  magrittr     * 2.0.1      2020-11-17 [1]
##  MASS           7.3-51.5   2019-12-20 [3]
##  Matrix         1.2-18     2019-11-27 [3]
##  memoise        1.1.0      2017-04-21 [2]
##  metap          1.2        2019-12-08 [2]
##  mnormt         1.5-5      2016-10-15 [2]
##  modelr         0.1.8      2020-05-19 [1]
##  multcomp       1.4-12     2020-01-10 [2]
##  multtest       2.42.0     2019-10-29 [2]
##  munsell        0.5.0      2018-06-12 [2]
##  mutoss         0.1-12     2017-12-04 [2]
##  mvtnorm        1.0-12     2020-01-09 [2]
##  nlme           3.1-143    2019-12-10 [3]
##  npsurv         0.4-0      2017-10-14 [2]
##  numDeriv       2016.8-1.1 2019-06-06 [2]
##  pbapply        1.4-2      2019-08-31 [2]
##  pillar         1.4.6      2020-07-10 [1]
##  pkgbuild       1.0.6      2019-10-09 [2]
##  pkgconfig      2.0.3      2019-09-22 [1]
##  pkgload        1.0.2      2018-10-29 [2]
##  plotly         4.9.1      2019-11-07 [2]
##  plotrix        3.7-7      2019-12-05 [2]
##  plyr           1.8.5      2019-12-10 [2]
##  png            0.1-7      2013-12-03 [2]
##  prettyunits    1.1.1      2020-01-24 [1]
##  processx       3.4.2      2020-02-09 [1]
##  progeny      * 1.11.3     2020-10-22 [1]
##  ps             1.3.2      2020-02-13 [1]
##  purrr        * 0.3.4      2020-04-17 [1]
##  R.methodsS3    1.7.1      2016-02-16 [2]
##  R.oo           1.23.0     2019-11-03 [2]
##  R.utils        2.9.2      2019-12-08 [2]
##  R6             2.4.1      2019-11-12 [1]
##  RANN           2.6.1      2019-01-08 [2]
##  rappdirs       0.3.1      2016-03-28 [2]
##  RColorBrewer   1.1-2      2014-12-07 [2]
##  Rcpp           1.0.4      2020-03-17 [1]
##  RcppAnnoy      0.0.16     2020-03-08 [1]
##  RcppParallel   4.4.4      2019-09-27 [2]
##  Rdpack         0.11-1     2019-12-14 [2]
##  readr        * 1.4.0      2020-10-05 [1]
##  readxl       * 1.3.1      2019-03-13 [2]
##  rematch        1.0.1      2016-04-21 [2]
##  remotes        2.1.0      2019-06-24 [2]
##  reprex         0.3.0      2019-05-16 [2]
##  reshape2       1.4.3      2017-12-11 [2]
##  reticulate     1.14       2019-12-17 [2]
##  rlang          0.4.8      2020-10-08 [1]
##  rmarkdown      2.0        2019-12-12 [2]
##  ROCR           1.0-7      2015-03-26 [2]
##  rprojroot      1.3-2      2018-01-03 [2]
##  rstudioapi     0.11       2020-02-07 [1]
##  rsvd           1.0.3      2020-02-17 [1]
##  Rtsne          0.15       2018-11-10 [2]
##  rvest          0.3.6      2020-07-25 [1]
##  sandwich       2.5-1      2019-04-06 [2]
##  scales         1.1.0      2019-11-18 [2]
##  sctransform    0.2.1      2019-12-17 [2]
##  SDMTools       1.1-221.2  2019-11-30 [2]
##  sessioninfo    1.1.1      2018-11-05 [2]
##  Seurat       * 3.1.2      2019-12-12 [2]
##  sn             1.5-4      2019-05-14 [2]
##  stringi        1.5.3      2020-09-09 [1]
##  stringr      * 1.4.0      2019-02-10 [1]
##  survival       3.1-8      2019-12-03 [3]
##  testthat       2.3.2      2020-03-02 [1]
##  TFisher        0.2.0      2018-03-21 [2]
##  TH.data        1.0-10     2019-01-21 [2]
##  tibble       * 3.0.4      2020-10-12 [1]
##  tidyr        * 1.1.2      2020-08-27 [1]
##  tidyselect     1.1.0      2020-05-11 [1]
##  tidyverse    * 1.3.0      2019-11-21 [2]
##  tsne           0.1-3      2016-07-15 [2]
##  usethis        1.5.1      2019-07-04 [2]
##  uwot           0.1.5      2019-12-04 [2]
##  vctrs          0.3.5      2020-11-17 [1]
##  viridis      * 0.5.1      2018-03-29 [2]
##  viridisLite  * 0.3.0      2018-02-01 [2]
##  withr          2.3.0      2020-09-22 [1]
##  xfun           0.12       2020-01-13 [2]
##  xml2           1.3.2      2020-04-23 [1]
##  yaml           2.2.1      2020-02-01 [1]
##  zoo            1.8-7      2020-01-10 [2]
##  source                                 
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.3)                         
##  CRAN (R 3.6.2)                         
##  Bioconductor                           
##  Bioconductor                           
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  Github (clauswilke/colorblindr@1ac3d4d)
##  R-Forge (R 3.6.2)                      
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.3)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.3)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.3)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.3)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  Bioconductor                           
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.3)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  Github (saezlab/progeny@94bea1c)       
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.3)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.3)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.3)                         
##  CRAN (R 3.6.3)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.3)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
## 
## [1] /home/uhlitzf/R/lib
## [2] /usr/local/lib/R/site-library
## [3] /usr/local/lib/R/library